home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / clang / jcool01.zip / MATRIX.C < prev    next >
C/C++ Source or Header  |  1992-09-25  |  26KB  |  634 lines
  1. //
  2. // Copyright (C) 1991 Texas Instruments Incorporated.
  3. // Copyright (C) 1992 General Electric Company.
  4. //
  5. // Permission is granted to any individual or institution to use, copy, modify,
  6. // and distribute this software, provided that this complete copyright and
  7. // permission notice is maintained, intact, in all copies and supporting
  8. // documentation.
  9. //
  10. // Texas Instruments Incorporated, General Electric Company,
  11. // provides this software "as is" without express or implied warranty.
  12. //
  13.  
  14. #include <cool/Matrix.h>            // header of matrix and its envelope
  15.  
  16. //## hack to workaround BC++ 3.1 Envelope bug
  17. #define CoolEnvelope CoolEnvelope_Matrix
  18.  
  19. // compare_s -- Pointer operator== function
  20. template<class Type> 
  21. Boolean (*CoolMatrix<Type>::compare_s)(const Type&, const Type&);
  22.  
  23.  
  24. // Matrix -- constructor specifiying size of Matrix
  25. // Input:    Row, Column parameters
  26. // Output:   None
  27.  
  28. template<class Type> 
  29. CoolMatrix<Type>::CoolMatrix(unsigned int rows, unsigned int cols)
  30. : CoolBase_Matrix(rows, cols) {
  31.   this->data = new Type*[rows];            // Allocate the row memory
  32.   Type* columns = new Type[cols*rows];        // Allocate the array of elmts
  33.   for (int i = 0; i < rows; i ++)        // For each row in the Matrix
  34.     this->data[i] = &columns[i*cols];        // Fill in address of row
  35.   if (this->compare_s == NULL)            // If not set yet
  36.     this->compare_s = &CoolMatrix_is_data_equal; // Default is_equal
  37. }
  38.  
  39.  
  40. // Matrix -- constructor specifiying size of matrix and initial value
  41. // Input:    Row, Column parameters and initial value
  42. // Output:   None
  43.  
  44. template<class Type> 
  45. CoolMatrix<Type>::CoolMatrix(unsigned int rows, unsigned int cols,
  46.                     const Type& value)
  47. : CoolBase_Matrix(rows, cols) {
  48.   this->data = new Type*[rows];            // Allocate the row memory
  49.   Type* columns = new Type[cols*rows];        // Allocate the array of elmts
  50.   for (int i = 0; i < rows; i ++) {        // For each row in the Matrix
  51.     this->data[i] = &columns[i*cols];        // Fill in address of row
  52.     for (int j = 0; j < cols; j++)        // For each element in column
  53.       this->data[i][j] = value;            // Assign initial value
  54.   }
  55.   if (this->compare_s == NULL)            // If not set yet
  56.     this->compare_s = &CoolMatrix_is_data_equal; // Default is_equal
  57. }
  58.  
  59. // Matrix -- constructor specifiying size of matrix and initial values
  60. // Input:    Rows, Column parameters and initial values
  61. // Output:   None
  62. // Note: Arguments in ... can only be pointers, primitive types like int, double,
  63. //       and NOT OBJECTS, passed by reference or value, like vectors, matrices;
  64. //       because constructors must be known and called at compile time!!!
  65. //       Cannot have char in ..., because char is 1 instead of 4 bytes, and 
  66. //       va_arg expects sizeof(Type) a multiple of 4 bytes.
  67.  
  68. template<class Type> 
  69. CoolMatrix<Type>::CoolMatrix(unsigned int rows, unsigned int cols, int n,
  70.                     Type v00 ...)
  71. : CoolBase_Matrix(rows, cols) {
  72. #if ERROR_CHECKING
  73.   if (((sizeof(Type) % 4) != 0) ||        // Cause alignment problems
  74.       (sizeof(Type) > 8))            // User defined classes?
  75.     this->va_arg_error(#Type, sizeof(Type));    // So, cannot use this constructor
  76. #endif
  77.   this->data = new Type*[rows];            // Allocate the row memory
  78.   Type* columns = new Type[cols*rows];        // Allocate the array of elmts
  79.   for (int i = 0; i < rows; i ++)        // For each row in the Matrix
  80.     this->data[i] = &columns[i*cols];        // Fill in address of row
  81.   if (n > 0) {                    // If user specified values
  82.     va_list argp;                // Declare argument list
  83.     va_start (argp, v00);                   // Initialize macro
  84.     for (i = 0; i < rows && n; i++)        // For remaining values given
  85.       for (int j = 0; j < cols && n; j++, n--)    // Moving sequentially in Matrix
  86.     if ((i == 0) && (j == 0))
  87.       this->data[0][0] = v00;               // Hack for v00 ...
  88.     else
  89.       this->data[i][j] = va_arg(argp, Type); // Extract and assign
  90.     va_end(argp);
  91.   }
  92.   if (this->compare_s == NULL)            // If not set yet
  93.     this->compare_s = &CoolMatrix_is_data_equal; // Default is_equal
  94. }
  95.  
  96. // Matrix -- Construct a matrix from a block array of data, stored row-wise.
  97. // Input     number of rows and columns, and array of r*c data.
  98. // Ouput     None
  99.  
  100. template<class Type>
  101. CoolMatrix<Type>::CoolMatrix(unsigned int rows, unsigned int cols,
  102.                     const Type* data_block) 
  103. : CoolBase_Matrix(rows,cols) {
  104.   this->data = new Type*[num_rows];        // Allocate the row memory
  105.   int n = num_rows * num_cols;            // Number of elements 
  106.   Type* columns = new Type[n];            // Allocate the array of elmts
  107.   for (int d = 0; d < n; d++)            // Copy all the data elements
  108.     columns[d] = data_block[d];
  109.   for (int i = 0; i < num_rows; i ++)         // For each row in the Matrix
  110.     this->data[i] = &columns[i*num_cols];    // Fill in address of row
  111. }
  112.  
  113.  
  114. // Matrix -- Copy constructor
  115. // Input:    Matrix reference
  116. // Output:   None
  117.  
  118. template<class Type> 
  119. CoolMatrix<Type>::CoolMatrix(const CoolMatrix<Type>& m)
  120. : CoolBase_Matrix(m) {
  121.   this->data = new Type*[this->num_rows];    // Allocate the row memory
  122.   Type* columns = new Type[num_cols*num_rows];    // Allocate the array of elmts
  123.   for (int i = 0; i < num_rows; i ++) {        // For each row in the Matrix
  124.     this->data[i] = &columns[i*num_cols];    // Fill in address of row
  125.     for (int j = 0; j < this->num_cols; j++)    // For each element in column
  126.       this->data[i][j] = m.data[i][j];        // Copy value
  127.   }
  128. }
  129.  
  130.  
  131. // ~Matrix -- Destructor for Matrix class that frees up storage
  132. // Input:     *this
  133. // Output:    None
  134.  
  135. template<class Type> 
  136. CoolMatrix<Type>::~CoolMatrix() {
  137.   delete [] this->data[0];            // Free up the array of elmts
  138.   delete [] this->data;                // Free up the row memory
  139. }
  140.  
  141. // fill -- Set all elements of a matrix to a specified fill value
  142. // Input:  this*, reference to fill value
  143. // Output: None
  144.  
  145. template<class Type> 
  146. void CoolMatrix<Type>::fill (const Type& value) {
  147.   for (int i = 0; i < this->num_rows; i++)    // For each row in the Matrix
  148.     for (int j = 0; j < this->num_cols; j++)    // For each element in column
  149.       this->data[i][j] = value;            // Assign fill value
  150. }
  151.  
  152. // operator= -- Overload the assignment operator to assign a single 
  153. //              value to the elements of a Matrix. 
  154. // Input:       *this, reference to a value
  155. // Output:      Reference to updated Matrix object
  156.  
  157. template<class Type> 
  158. CoolMatrix<Type>& CoolMatrix<Type>::operator= (const Type& value) {
  159.   for (int i = 0; i < this->num_rows; i++)    // For each row in Matrix
  160.     for (int j = 0; j < this->num_cols; j++)    // For each column in Matrix
  161.       this->data[i][j] = value;            // Assign value
  162.   return *this;                    // Return Matrix reference
  163. }
  164.  
  165.  
  166. // operator= -- Overload the assignment operator to copy the elements
  167. //              in one Matrix to another. The existing storage for the 
  168. //              destination matrix is freed up and new storage of the same 
  169. //              size as the source is allocated.
  170. // Input:       *this, reference to Matrix
  171. // Output:      Reference to copied Matrix object
  172.  
  173. template<class Type> 
  174. CoolMatrix<Type>& CoolMatrix<Type>::operator= (const CoolMatrix<Type>& m) {
  175.   if (this != &m) {                // make sure *this != m
  176.     if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) {
  177.       delete [] this->data[0];            // Free up the array of elmts
  178.       delete [] this->data;            // Free up the row memory
  179.       this->num_rows = m.num_rows;        // Copy rows
  180.       this->num_cols = m.num_cols;        // Copy columns 
  181.       this->data = new Type*[this->num_rows];    // Allocate the rows
  182.       Type* columns = new Type[num_cols*num_rows]; // Allocate the columns
  183.       for (int i = 0; i < this->num_rows; i++)       // For each row
  184.     this->data[i] = &columns[i*num_cols];       // Fill in address of row
  185.     }
  186.     for (int i = 0; i < this->num_rows; i++)    // For each row in the Matrix
  187.       for (int j = 0; j < this->num_cols; j++)    // For each element in column
  188.     this->data[i][j] = m.data[i][j];    // Copy value
  189.   }
  190.   return *this;                    // Return Matrix reference
  191. }
  192.  
  193. // operator== -- Compare the elements of two Matrices of Type Type using
  194. //               the Compare pointer to funtion (default is ==). If one 
  195. //               Matrix has more rows and/or columns than the other, the
  196. //               result is FALSE
  197. // Input:        Reference to Matrix of Type Type
  198. // Output:       TRUE/FALSE
  199.  
  200. template<class Type> 
  201. Boolean CoolMatrix<Type>::operator== (const CoolMatrix<Type>& m) const {
  202.   if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
  203.     return FALSE;                            // Then not equal
  204.   for (int i = 0; i < this->num_rows; i++)                // For each row
  205.     for (int j = 0; j < this->num_cols; j++)                // For each columne
  206.       if ((*this->compare_s)(this->data[i][j],m.data[i][j]) == FALSE) // Same?
  207.     return FALSE;                              // Then no match
  208.   return TRUE;                              // Else same, so return TRUE
  209. }
  210.  
  211.  
  212. // is_data_equal -- Default data comparison function if user has not provided
  213. //                  another one. Note that this is not inline because we need
  214.     //                  to take the address of it for the compare static variable
  215. // Input:           Two Type references
  216. // Output:          TRUE/FALSE
  217.  
  218. template<class Type>
  219.   Boolean CoolMatrix_is_data_equal (const Type& t1, const Type& t2) {
  220.     return (t1 == t2);
  221.   }
  222.  
  223.  
  224. // set_compare -- Specify the comparison function to be used
  225. //                in logical tests of vector elements
  226. // Input:         Pointer to a compare function
  227. // Output:        None
  228.  
  229. template<class Type> 
  230. void CoolMatrix<Type>::set_compare (/*Compare##*/Boolean (*c)(const Type&, const Type&)) {
  231.   if (c == NULL)                // If no argument supplied
  232.     this->compare_s = &CoolMatrix_is_data_equal; // Default is_equal
  233.   else
  234.     this->compare_s = c;            // Else set to user function
  235. }
  236.  
  237. // operator<< -- Overload the output operator to print a matrix
  238. // Input:        ostream reference, Matrix reference
  239. // Output:       ostream reference
  240.  
  241. template<class Type>
  242.   ostream& operator<< (ostream& s, const CoolMatrix<Type>& m) {
  243.     for (int i = 0; i < m.rows(); i++) {    // For each row in matrix
  244.       for (int j = 0; j < m.columns(); j++)    // For each column in matrix
  245.     s << m.data[i][j] << " ";        // Output data element
  246.       s << "\n";                    // Output newline
  247.     }
  248.     return (s);                    // Return ostream reference
  249. }
  250.  
  251.  
  252.  
  253. // operator+= -- Destructive matrix addition of a scalar.
  254. // Input:        this*, scalar value
  255. // Output:       New matrix reference
  256.  
  257. template<class Type> 
  258. CoolMatrix<Type>& CoolMatrix<Type>::operator+= (const Type& value) {
  259.   for (int i = 0; i < this->num_rows; i++)    // For each row
  260.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  261.       this->data[i][j] += value;        // Add scalar
  262.   return *this;
  263. }
  264.  
  265. // operator*= -- Destructive matrix multiplication by a scalar.
  266. // Input:        this*, scalar value
  267. // Output:       New matrix reference
  268.  
  269.  
  270. template<class Type> 
  271. CoolMatrix<Type>& CoolMatrix<Type>::operator*= (const Type& value) {
  272.   for (int i = 0; i < this->num_rows; i++)    // For each row
  273.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  274.       this->data[i][j] *= value;        // Multiply by scalar
  275.   return *this;
  276. }
  277.  
  278. // operator/= -- Destructive matrix division by a scalar.
  279. // Input:        this*, scalar value
  280. // Output:       New matrix reference
  281.  
  282. template<class Type> 
  283. CoolMatrix<Type>& CoolMatrix<Type>::operator/= (const Type& value) {
  284.   for (int i = 0; i < this->num_rows; i++)    // For each row
  285.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  286.       this->data[i][j] /= value;        // division by scalar
  287.   return *this;
  288. }
  289.  
  290.  
  291. // operator+= -- Destructive matrix addition with assignment. Note that the
  292. //               dimensions of each matrix must be identical
  293. // Input:        this*, matrix reference
  294. // Output:       Updated this* matrix reference
  295.  
  296. template<class Type> 
  297. CoolMatrix<Type>& CoolMatrix<Type>::operator+= (const CoolMatrix<Type>& m) {
  298.   if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
  299.     this->dimension_error ("operator+=", "Type", 
  300.                this->num_rows, this->num_cols, m.num_rows, m.num_cols);
  301.   for (int i = 0; i < this->num_rows; i++)    // For each row
  302.     for (int j = 0; j < this->num_cols; j++)    // For each element in column
  303.       this->data[i][j] += m.data[i][j];        // Add elements
  304.   return *this;
  305. }
  306.  
  307.  
  308. // operator-= -- Destructive matrix subtraction with assignment. Note that the
  309. //               dimensions of each matrix must be identical
  310. // Input:        this*, matrix reference
  311. // Output:       Updated this* matrix reference
  312.  
  313. template<class Type> 
  314. CoolMatrix<Type>& CoolMatrix<Type>::operator-= (const CoolMatrix<Type>& m) {
  315.   if (this->num_rows != m.num_rows || this->num_cols != m.num_cols) // Size?
  316.     this->dimension_error ("operator-=", "Type", 
  317.                this->num_rows, this->num_cols, m.num_rows, m.num_cols);
  318.   int i, j;
  319.   for (i = 0; i < this->num_rows; i++)
  320.     for (j = 0; j < this->num_cols; j++)
  321.       this->data[i][j] -= m.data[i][j];
  322.   return *this;
  323. }
  324.  
  325. // operator* -- Non Destructive matrix multiply 
  326. //               num_cols of first matrix must match num_rows of second matrix.
  327. // Input:        two matrix references
  328. // Output:       New matrix containing the product.
  329.  
  330. template<class Type>
  331. CoolEnvelope< CoolMatrix<Type> > operator* (const CoolMatrix<Type>& m1,
  332.                  const CoolMatrix<Type>& m2) {
  333.   if (m1.num_cols != m2.num_rows)        // dimensions do not match?
  334.     m1.dimension_error ("operator*=", "Type", 
  335.             m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  336.   CoolMatrix<Type> temp(m1.num_rows, m2.num_cols); // Temporary to store product
  337.   for (int i = 0; i < m1.num_rows; i++) {    // For each row
  338.     for (int j = 0; j < m2.num_cols; j++) {    // For each element in column
  339.       Type sum = 0;
  340.       for (int k = 0; k < m1.num_cols; k++)    // Loop over column values
  341.     sum += (m1.data[i][k] * m2.data[k][j]); // Multiply
  342.       temp(i,j) = sum;
  343.     }
  344.   }
  345.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  346.   return result;                     // copy of envelope
  347. }
  348.  
  349. // operator- -- Non-destructive matrix negation. a = -b;
  350. // Input:       this*
  351. // Output:      New matrix 
  352.  
  353. template<class Type> 
  354. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator- () const {
  355.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  356.   int i, j;
  357.   for (i = 0; i < this->num_rows; i++)
  358.     for (j = 0; j < this->num_cols; j++)
  359.       temp.data[i][j] = - this->data[i][j];
  360.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  361.   return result;                     // copy of envelope
  362. }
  363.  
  364. // operator+ -- Non-destructive matrix addition of a scalar.
  365. // Input:       this*, scalar value
  366. // Output:      New matrix 
  367.  
  368. template<class Type> 
  369. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator+ (const Type& value) const {
  370.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  371.   for (int i = 0; i < this->num_rows; i++)    // For each row
  372.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  373.       temp.data[i][j] = (this->data[i][j] + value); // Add scalar
  374.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  375.   return result;                     // copy of envelope
  376. }
  377.  
  378.  
  379. // operator* -- Non-destructive matrix multiply by a scalar.
  380. // Input:       this*, scalar value
  381. // Output:      New matrix 
  382.  
  383. template<class Type> 
  384. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator* (const Type& value) const {
  385.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  386.   for (int i = 0; i < this->num_rows; i++)    // For each row
  387.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  388.       temp.data[i][j] = (this->data[i][j] * value); // Multiply
  389.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  390.   return result;                     // copy of envelope
  391. }
  392.  
  393.  
  394. // operator/ -- Non-destructive matrix division by a scalar.
  395. // Input:       this*, scalar value
  396. // Output:      New matrix 
  397.  
  398. template<class Type> 
  399. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::operator/ (const Type& value) const {
  400.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  401.   for (int i = 0; i < this->num_rows; i++)    // For each row
  402.     for (int j = 0; j < this->num_cols; j++)     // For each element in column
  403.       temp.data[i][j] = (this->data[i][j] / value); // Divide
  404.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  405.   return result;                     // copy of envelope
  406. }
  407.  
  408.  
  409. ////--------------------------- Additions ------------------------------------
  410.  
  411.  
  412. // transpose -- Return the transpose of this matrix.
  413. // Input:       this*
  414. // Ouput:       New matrix
  415.  
  416. template<class Type> 
  417. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::transpose() const {
  418.   CoolMatrix<Type> temp(this->num_cols, this->num_rows);
  419.   int i, j;
  420.   for (i = 0; i < this->num_cols; i++)
  421.     for (j = 0; j < this->num_rows; j++)
  422.       temp.data[i][j] = this->data[j][i];
  423.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  424.   return result;                     // copy of envelope
  425. }
  426.  
  427.  
  428. // abs -- Return the matrix of the absolute values.
  429. // Input:       this*
  430. // Ouput:       New matrix
  431.  
  432. template<class Type> 
  433. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::abs() const {
  434.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  435.   int i, j;
  436.   for (i = 0; i < this->num_rows; i++)
  437.     for (j = 0; j < this->num_cols; j++)
  438.       if (this->data[i][j] < 0)
  439.     temp.data[i][j] = - this->data[i][j];
  440.       else
  441.     temp.data[i][j] = this->data[i][j];
  442.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  443.   return result;                     // copy of envelope
  444. }
  445.  
  446. // sign -- Return the matrix whose elements are either -1,1 or 0
  447. // depending on whether the corresponding values are negative, positive, or 0.
  448. // Input:       this*
  449. // Ouput:       New matrix
  450.  
  451. template<class Type> 
  452. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::sign() const {
  453.   CoolMatrix<Type> temp(this->num_rows, this->num_cols);
  454.   int i, j;
  455.   for (i = 0; i < this->num_rows; i++)
  456.     for (j = 0; j < this->num_cols; j++)
  457.       if (this->data[i][j] == 0)            // test fuzz equality to 0
  458.     temp.data[i][j] = 0;            // first.
  459.       else
  460.     if (this->data[i][j] < 0)
  461.       temp.data[i][j] = -1;
  462.     else
  463.       temp.data[i][j] = 1;
  464.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  465.   return result;                     // copy of envelope
  466. }
  467.  
  468. // element_product -- return the matrix whose elements are the products 
  469. // Input:     2 matrices m1, m2 by reference
  470. // Output:    New matrix, whose elements are m1[ij]*m2[ij].
  471.  
  472. template<class Type>
  473. CoolEnvelope< CoolMatrix<Type> > element_product (const CoolMatrix<Type>& m1, const CoolMatrix<Type>& m2) {
  474.   if (m1.num_rows != m2.num_rows || m1.num_cols != m2.num_cols) // Size?
  475.     m1.dimension_error ("element_product", "Type", 
  476.             m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  477.   CoolMatrix<Type> temp(m1.num_rows, m1.num_cols);
  478.   int i, j;
  479.   for (i = 0; i < m1.num_rows; i++)
  480.     for (j = 0; j < m1.num_cols; j++)
  481.       temp.data[i][j] = m1.data[i][j] * m2.data[i][j];
  482.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  483.   return result;                     // copy of envelope
  484. }
  485.  
  486. // element_quotient -- return the matrix whose elements are the quotients 
  487. // Input:     2 matrices m1, m2 by reference
  488. // Output:    New matrix, whose elements are m1[ij]/m2[ij].
  489.  
  490. template<class Type>
  491. CoolEnvelope< CoolMatrix<Type> > element_quotient (const CoolMatrix<Type>& m1, const CoolMatrix<Type>& m2) {
  492.   if (m1.num_rows != m2.num_rows || m1.num_cols != m2.num_cols) // Size?
  493.     m1.dimension_error ("element_quotient", "Type", 
  494.             m1.num_rows, m1.num_cols, m2.num_rows, m2.num_cols);
  495.   CoolMatrix<Type> temp(m1.num_rows, m1.num_cols);
  496.   int i, j;
  497.   for (i = 0; i < m1.num_rows; i++)
  498.     for (j = 0; j < m1.num_cols; j++)
  499.       temp.data[i][j] = m1.data[i][j] / m2.data[i][j];
  500.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  501.   return result;                     // copy of envelope
  502. }
  503.  
  504. // update -- replace a submatrix of this, by the actual argument.
  505. // Input:       *this, starting corner specified by top and left.
  506. // Ouput:       mutated reference.
  507.  
  508. template<class Type> 
  509. CoolMatrix<Type>& CoolMatrix<Type>::update (const CoolMatrix<Type>& m, 
  510.                         unsigned int top, unsigned int left) {
  511.   unsigned int bottom = top + m.num_rows;
  512.   unsigned int right = left + m.num_cols;
  513.   if (this->num_rows < bottom || this->num_cols < right)
  514.     this->dimension_error ("update", "Type", 
  515.                bottom-top, right-left, m.num_rows, m.num_cols);
  516.   int i, j;
  517.   for (i = top; i < bottom; i++)
  518.     for (j = left; j < right; j++)
  519.       this->data[i][j] = m.data[i-top][j-left];
  520.   return *this;
  521. }
  522.  
  523.  
  524. // extract -- Return a submatrix specified by the top-left corner and size.
  525. // Input:       *this, starting corner specified by top and left, and size.
  526. // Ouput:       new matrix
  527.  
  528. template<class Type> 
  529. CoolEnvelope< CoolMatrix<Type> > CoolMatrix<Type>::extract (unsigned int rows, unsigned int cols, unsigned int top, unsigned int left) const{
  530.   unsigned int bottom = top + rows;
  531.   unsigned int right = left + cols;
  532.   if ((this->num_rows < bottom) || (this->num_cols < right))
  533.     this->dimension_error ("extract", "Type", 
  534.                bottom-top, right-left, rows, cols);
  535.   CoolMatrix<Type> temp(rows, cols);
  536.   for (int i = 0; i < rows; i++)        // actual copy of all elements
  537.     for (int j = 0; j < cols ; j++)        // in submatrix
  538.       temp.data[i][j] = data[top+i][left+j];
  539.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  540.   return result;                     // copy of envelope
  541. }
  542.  
  543. // determinant -- Determinant of a square matrix using Cramer's rule.
  544. //              Signal Error exception if the matrix is not square.
  545.  
  546. template<class Type>
  547. Type CoolMatrix<Type>::determinant () const {
  548.   if (this->num_rows != this->num_cols || this->num_rows < 2)
  549.     this->dimension_error ("determinant", "Type",
  550.                this->num_rows, this->num_cols, 
  551.                this->num_rows, this->num_cols);
  552.   int n = this->num_rows, r, i, j;
  553.   Type det = 0, prod;
  554.   if (n == 2) {
  555.     det = (this->data[0][0] * this->data[1][1] - // border case of 2x2 matrix
  556.        this->data[0][1] * this->data[1][0]);
  557.   } else {
  558.     for (r = 0; r < n; r++) {            // compute sum of products
  559.       prod = 1;                    // along diagonals
  560.       for (i = r, j = 0; i < n; i++, j++)    // top-left to bottom-right
  561.     prod *= this->data[i][j];
  562.       for (i = 0; j < n; i++, j++)
  563.     prod *= this->data[i][j];
  564.       det += prod;                // coeft = +1
  565.     }
  566.     int e = n-1;                // index of last row/col
  567.     for (r = 0; r < n; r++) {            // compute sum of products
  568.       prod = 1;                    // along diagonals
  569.       for (i = r, j = e; i < n; i++, j--)    // top-right to bottom-left
  570.     prod *= this->data[i][j];
  571.       for (i=0; j >= 0; i++, j--)
  572.     prod *= this->data[i][j];
  573.     det -= prod;                // coeft = -1
  574.     }
  575.   }
  576.   return det;
  577. }
  578.  
  579. // dot_product -- Return the dot product of the row or column vectors
  580. // Input:       2 vectors by reference
  581. // Ouput:       dot product value
  582.  
  583. template<class Type>
  584.   Type dot_product (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  585.     if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols) // Size?
  586.       v1.dimension_error ("dot_product", "Type", 
  587.               v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  588.     Type dot = 0;
  589.     int i, j;
  590.     for (i = 0; i < v1.num_rows; i++)
  591.       for (j = 0; j < v1.num_cols; j++)        // generalized dot-product
  592.     dot += v1.data[i][j] * v2.data[i][j];    // of matrices
  593.     return dot;
  594. }
  595.  
  596. // cross_2d -- Return the 2X1 cross-product of 2 2d-vectors
  597. // Input:       2 vectors by reference
  598. // Ouput:       cross product value
  599.  
  600. template<class Type>
  601. Type cross_2d (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  602.   if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols)
  603.     v1.dimension_error ("cross_2d", "Type", 
  604.             v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  605.   CoolMatrix<Type>& m1 = (CoolMatrix<Type>&) v1; // cast away const.
  606.   CoolMatrix<Type>& m2 = (CoolMatrix<Type>&) v2; 
  607.   return (m1.x() * m2.y()            // work for both col/row
  608.       -                    // representation.
  609.       m1.y() * m2.x());
  610. }
  611.  
  612. // cross_3d -- Return the 3X1 cross-product of 2 3d-vectors
  613. // Input:       2 vectors by reference
  614. // Ouput:       3d cross product vector
  615.  
  616. template<class Type>
  617. CoolEnvelope< CoolMatrix<Type> > cross_3d (const CoolMatrix<Type>& v1, const CoolMatrix<Type>& v2) {
  618.   if (v1.num_rows != v2.num_rows || v1.num_cols != v2.num_cols)
  619.     v1.dimension_error ("cross_3d", "Type", 
  620.             v1.num_rows, v1.num_cols, v2.num_rows, v2.num_cols);
  621.   CoolMatrix<Type> temp(v1.num_rows, v1.num_cols);
  622.   CoolMatrix<Type>& m1 = (CoolMatrix<Type>&) v1; // cast away const.
  623.   CoolMatrix<Type>& m2 = (CoolMatrix<Type>&) v2; 
  624.   temp.x() = m1.y() * m2.z() - m1.z() * m2.y(); // work for both col/row
  625.   temp.y() = m1.z() * m2.x() - m1.x() * m2.z(); // representation
  626.   temp.z() = m1.x() * m2.y() - m1.y() * m2.x();
  627.   CoolEnvelope< CoolMatrix<Type> >& result = (CoolEnvelope< CoolMatrix<Type> >&) temp; // same physical object
  628.   return result;                     // copy of envelope
  629. }
  630.  
  631. //## hack to workaround BC++ 3.1 Envelope bug
  632. #undef CoolEnvelope
  633.  
  634.